#include "MeshProcess.hpp"

void INN_IEN_2d_setup( const vector<double> sKnots,
    const vector<double> tKnots,
    const int sDegree,
    const int tDegree,
    int ** &INN,
    int ** &IEN,
    double * &hx,
    double * &hy,
    int &nElem,
    int &nElem_x,
    int &nElem_y,
    int &nFunc,
    int &nFunc_x,
    int &nFunc_y,
    int &nLocBas )
{
  // n and m are the number of function in xi and eta direction
  int n, m; 
  n = sKnots.size() - sDegree - 1;
  m = tKnots.size() - tDegree - 1;

  nFunc_x = n;
  nFunc_y = m;

  nFunc = n * m;
  nElem = (n-sDegree)*(m-tDegree);

  nElem_x = n - sDegree;
  nElem_y = m - tDegree;

  nLocBas = (sDegree+1)*(tDegree+1);

  IEN = new int * [nElem];
  int ** IENv;
  IENv = new int * [nElem];

  for(int i=0; i<nElem; i++)
  {
    IEN[i] = new int[nLocBas];
    IENv[i] = new int[nLocBas];
  } 

  INN = new int * [nFunc];
  for(int i=0; i<nFunc; i++)
  {
    INN[i] = new int[2];
  }

  hx = new double [nElem_x];
  hy = new double [nElem_y];

  for(int i=0; i<nElem_x; i++)
  {
    hx[i] = sKnots[sDegree+i+1] - sKnots[sDegree+i];
  }

  for(int i=0; i<nElem_y; i++)
  {
    hy[i] = tKnots[tDegree+i+1] - tKnots[tDegree+i];
  }


  int A = -1;
  int e = -1;
  for(int j = 0; j < m; j++)
  {
    for(int i = 0; i < n; i++)
    {
      A = A + 1;
      INN[A][0] = i;
      INN[A][1] = j;
      if((i>=sDegree) && (j>=tDegree))
      {
        e = e+1;
        for(int jloc = 0; jloc<=tDegree; jloc++)
        {
          for(int iloc = 0; iloc<=sDegree; iloc++)
          {
            int B = A - jloc*n - iloc;
            int b = jloc * (sDegree + 1) + iloc;
            IENv[e][b] = B;
          }
        }
      }
    }
  }


  for(int i=0; i<nElem; i++)
  {
    for(int j=0; j<nLocBas; j++)
    {
      IEN[i][j] = IENv[i][nLocBas-1-j];
    }
  }

  for(int i=0; i<nElem; i++)
    delete [] IENv[i];
  delete [] IENv;

};

void INN_IEN_2d_cleanup( int ** &IEN, int ** &INN, 
                         int nElem, int nFunc,
                         double * &hx, double * &hy )
{
  for(int i=0; i<nElem; i++)
  {
    delete [] IEN[i];
  }
  for(int i=0; i<nFunc; i++)
    delete [] INN[i];

  delete [] IEN;
  delete [] INN;
  delete [] hx;
  delete [] hy;
};

void QuadratureRule(int nqp, double *&qp, double *&wq)
{
  qp = new double [nqp];
  wq = new double [nqp];
  switch (nqp)
  {
    case 1:
      qp[0] = 0.5;
      wq[0] = 1.0;
      break;
    case 2:
      qp[0] = 0.788675134594813;
      wq[0] = 0.5;
      qp[1] = 0.211324865405187;
      wq[1] = 0.5;
      break;
    case 3:
      qp[0] = 0.887298334620742;
      qp[1] = 0.5;
      qp[2] = 0.112701665379258;
      wq[0] = 0.277777777777777;
      wq[1] = 0.444444444444444;
      wq[2] = 0.277777777777777;
      break;
    case 4:
      qp[0] = 0.930568155797026;
      qp[1] = 0.669990521792428;
      qp[2] = 0.330009478207572;
      qp[3] = 0.069431844202974;

      wq[0] = 0.173927422568727;
      wq[1] = 0.326072577431273;
      wq[2] = 0.326072577431273;
      wq[3] = 0.173927422568727;
      break;
    default:
      cerr<<"ERROR: We have quadrature rule with up to 4 quadrature points"<<endl;
      assert(0);
  }
};

void QuadratureClean(double * &qp, double * &wq)
{
  delete [] qp;
  delete [] wq;
};






























